%% process_showbasis
function process_showbasis(basis_type, derivative_type, varargin)
[x, y] = meshgrid(0:0.0125:1, 0:0.0125:1);
z = basis_ref(x(:), y(:), basis_type, derivative_type);
if isempty(varargin)
    f_begin = 1;
    f_end = size(z,2);
else
    f_begin = varargin{1};
    f_end = varargin{1};
    if f_end > size(z,2)
        f_end = size(z,2);
    end
end
triangle_region = @(x,y) (x >= 0) & (x <= 1) & (y >= 0) & (y <= 1) & (y <= -x+1);
figure("WindowStyle", "docked");
for j = f_begin:f_end
    Z = reshape(z(:,j), size(x,1), size(x,2));
    Z(~triangle_region(x, y)) = NaN;
    surf(x, y, Z);
    hold on
end
xlabel("x");
ylabel("y");
axis("tight", "equal");
view(-45,60);
end

%% basis_ref
function N = basis_ref(x, y, basis_type, derivative_type)
if strcmp(basis_type, "P1")
    N = zeros(size(x,1), 3);
elseif strcmp(basis_type, "P1b")
    N = zeros(size(x,1), 4);
elseif strcmp(basis_type, "P2")
    N = zeros(size(x,1), 6);
elseif strcmp(basis_type, "P3")
    N = zeros(size(x,1), 10);
elseif strcmp(basis_type, "P4")
    N = zeros(size(x,1), 15);
elseif strcmp(basis_type, "P5")
    N = zeros(size(x,1), 21);
else
    error("Invalid basis type.");
end
if strcmp(basis_type, "P1")
    if strcmp(derivative_type, "x")
        N(:,1) = 1-x-y;
        N(:,2) = x;
        N(:,3) = y;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
    end
elseif strcmp(basis_type, "P1b")
    if strcmp(derivative_type, "x")
        N(:,1) = 1-x-y;
        N(:,2) = x;
        N(:,3) = y;
        N(:,4) = 27*x.*y.*N(:,1);
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
        N(:,4) = -27*y.*(2*x+y-1);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
        N(:,4) = -27*x.*(x+2*y-1);
    end
elseif strcmp(basis_type, "P2")
    lambda = 1-x-y;
    if strcmp(derivative_type, "x")
        N(:,1) = lambda.*(2*lambda-1);
        N(:,4) = 4*x.*lambda;
        N(:,2) = x.*(2*x-1);
        N(:,5) = 4*x.*y;
        N(:,3) = y.*(2*y-1);
        N(:,6) = 4*y.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = 1-4*lambda;
        N(:,4) = 4*(lambda-x);
        N(:,2) = 4*x-1;
        N(:,5) = 4*y;
        N(:,3) = 0;
        N(:,6) = -4*y;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = 1-4*lambda;
        N(:,4) = -4*x;
        N(:,2) = 0;
        N(:,5) = 4*x;
        N(:,3) = 4*y-1;
        N(:,6) = 4*(lambda-y);
    end
elseif strcmp(basis_type, "P3")
    lambda = 1-x-y;
    if strcmp(derivative_type, "x")
        N(:,1) = 0.5*lambda.*(-1+3*lambda).*(-2+3*lambda);
        N(:,4) = 4.5*lambda.*x.*(-1+3*lambda);
        N(:,7) = 4.5*lambda.*x.*(-1+3*x);
        N(:,2) = 0.5*x.*(-1+3*x).*(-2+3*x);
        N(:,5) = 4.5*x.*y.*(-1+3*x);
        N(:,8) = 4.5*x.*y.*(-1+3*y);
        N(:,3) = 0.5*y.*(-1+3*y).*(-2+3*y);
        N(:,6) = 4.5*lambda.*y.*(-1+3*y);
        N(:,9) = 4.5*lambda.*y.*(-1+3*lambda);
        N(:,10) = 27*x.*y.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1+9*lambda-13.5*lambda.^2;
        N(:,4) = 4.5*lambda.*(-1+3*lambda-6*x)+4.5*x;
        N(:,7) = 4.5*x.*(1+6*lambda-3*x)-4.5*lambda;
        N(:,2) = 1-9*x+13.5*x.^2;
        N(:,5) = 4.5*y.*(-1+6*x);
        N(:,8) = 4.5*y.*(-1+3*y);
        N(:,3) = 0;
        N(:,6) = -4.5*y.*(-1+3*y);
        N(:,9) = -4.5*y.*(-1+6*lambda);
        N(:,10) = 27*y.*(lambda-x);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1+9*lambda-13.5*lambda.^2;
        N(:,4) = -4.5*x.*(-1+6*lambda);
        N(:,7) = -4.5*x.*(-1+3*x);
        N(:,2) = 0;
        N(:,5) = 4.5*x.*(-1+3*x);
        N(:,8) = 4.5*x.*(-1+6*y);
        N(:,3) = 1-9*y+13.5*y.^2;
        N(:,6) = 4.5*y.*(1+6*lambda-3*y)-4.5*lambda;
        N(:,9) = 4.5*lambda.*(-1+3*lambda-6*y)+4.5*y;
        N(:,10) = 27*x.*(lambda-y);
    end
elseif strcmp(basis_type, "P4")
    lambda = 1-x-y;
    if strcmp(derivative_type, "x")
        N(:,1) = lambda.*(2*lambda-1).*(4*lambda-3).*(4*lambda-1)/3;
        N(:,4) = 16*x.*(4*lambda-1).*(2*lambda-1).*lambda/3;
        N(:,7) = 4*x.*(4*lambda-1).*(4*x-1).*lambda;
        N(:,10) = 16*x.*(4*x-1).*(2*x-1).*lambda/3;
        N(:,2) = x.*(2*x-1).*(4*x-3).*(4*x-1)/3;
        N(:,5) = 16*x.*y.*(4*x-1).*(2*x-1)/3;
        N(:,8) = 4*x.*y.*(4*x-1).*(4*y-1);
        N(:,11) = 16*x.*y.*(4*y-1).*(2*y-1)/3;
        N(:,3) = y.*(2*y-1).*(4*y-3).*(4*y-1)/3;
        N(:,6) = 16*lambda.*y.*(4*y-1).*(2*y-1)/3;
        N(:,9) = 4*lambda.*y.*(4*lambda-1).*(4*y-1);
        N(:,12) = (16*lambda.*y.*(4*lambda-1).*(2*lambda-1))/3;
        N(:,13) = 32*x.*y.*(4*lambda-1).*lambda;
        N(:,14) = 32*x.*y.*(4*x-1).*lambda;
        N(:,15) = 32*x.*y.*(4*y-1).*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -(2*lambda-1).*(4*lambda-3).*((4*lambda)/3-1/3)-(4*(2*lambda-1).*(4*lambda-3).*lambda)/3-4*(2*lambda-1).*((4*lambda)/3-1/3).*lambda-2*(4*lambda-3).*((4*lambda)/3-1/3).*lambda;
        N(:,12) = -(16*y.*(2*lambda-1/2).*(4*lambda-2))/3-(64*y.*(2*lambda-1/2).*lambda)/3-(32*y.*(4*lambda-2).*lambda)/3;
        N(:,9) = -4*y.*(4*lambda-1).*(4*y-1)-16*y.*(4*y-1).*lambda;
        N(:,6) = -(16*y.*(2*y-1/2).*(4*y-2))/3;
        N(:,3) = 0;
        N(:,4) = (16*(2*lambda-1/2).*(4*lambda-2).*lambda)/3-(64*x.*(2*lambda-1/2).*lambda)/3-(32*x.*(4*lambda-2).*lambda)/3-(16*x.*(2*lambda-1/2).*(4*lambda-2))/3;
        N(:,13) = 32*y.*(4*lambda-1).*lambda-32*x.*y.*(4*lambda-1)-128*x.*y.*lambda;
        N(:,15) = 32*y.*(4*y-1).*lambda-32*x.*y.*(4*y-1);
        N(:,11) = (16*y.*(2*y-1/2).*(4*y-2))/3;
        N(:,7) = 4*(4*lambda-1).*(4*x-1).*lambda-4*x.*(4*lambda-1).*(4*x-1)+16*x.*(4*lambda-1).*lambda-16*x.*(4*x-1).*lambda;
        N(:,14) = 32*y.*(4*x-1).*lambda-32*x.*y.*(4*x-1)+128*x.*y.*lambda;
        N(:,8) = 4*y.*(4*x-1).*(4*y-1)+16*x.*y.*(4*y-1);
        N(:,10) = (16*(2*x-1/2).*(4*x-2).*lambda)/3-(16*x.*(2*x-1/2).*(4*x-2))/3+(64*x.*(2*x-1/2).*lambda)/3+(32*x.*(4*x-2).*lambda)/3;
        N(:,5) = (16*y.*(2*x-1/2).*(4*x-2))/3+(64*x.*y.*(2*x-1/2))/3+(32*x.*y.*(4*x-2))/3;
        N(:,2) = (2*x-1).*(4*x-3).*((4*x)/3-1/3)+(4*x.*(2*x-1).*(4*x-3))/3+4*x.*(2*x-1).*((4*x)/3-1/3)+2*x.*(4*x-3).*((4*x)/3-1/3);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -(2*lambda-1).*(4*lambda-3).*((4*lambda)/3-1/3)-(4*(2*lambda-1).*(4*lambda-3).*lambda)/3-4*(2*lambda-1).*((4*lambda)/3-1/3).*lambda-2*(4*lambda-3).*((4*lambda)/3-1/3).*lambda;
        N(:,12) = (16*(2*lambda-1/2).*(4*lambda-2).*lambda)/3-(64*y.*(2*lambda-1/2).*lambda)/3-(32*y.*(4*lambda-2).*lambda)/3-(16*y.*(2*lambda-1/2).*(4*lambda-2))/3;
        N(:,9) = 4*(4*lambda-1).*(4*y-1).*lambda-4*y.*(4*lambda-1).*(4*y-1)+16*y.*(4*lambda-1).*lambda-16*y.*(4*y-1).*lambda;
        N(:,6) = (16*(2*y-1/2).*(4*y-2).*lambda)/3-(16*y.*(2*y-1/2).*(4*y-2))/3+(64*y.*(2*y-1/2).*lambda)/3+(32*y.*(4*y-2).*lambda)/3;
        N(:,3) = (2*y-1).*(4*y-3).*((4*y)/3-1/3)+(4*y.*(2*y-1).*(4*y-3))/3+4*y.*(2*y-1).*((4*y)/3-1/3)+2*y.*(4*y-3).*((4*y)/3-1/3);
        N(:,4) = -(16*x.*(2*lambda-1/2).*(4*lambda-2))/3-(64*x.*(2*lambda-1/2).*lambda)/3-(32*x.*(4*lambda-2).*lambda)/3;
        N(:,13) = 32*x.*(4*lambda-1).*lambda-32*x.*y.*(4*lambda-1)-128*x.*y.*lambda;
        N(:,15) = 32*x.*(4*y-1).*lambda-32*x.*y.*(4*y-1)+128*x.*y.*lambda;
        N(:,11) = (16*x.*(2*y-1/2).*(4*y-2))/3+(64*x.*y.*(2*y-1/2))/3+(32*x.*y.*(4*y-2))/3;
        N(:,7) = -4*x.*(4*lambda-1).*(4*x-1)-16*x.*(4*x-1).*lambda;
        N(:,14) = 32*x.*(4*x-1).*lambda-32*x.*y.*(4*x-1);
        N(:,8) = 4*x.*(4*x-1).*(4*y-1)+16*x.*y.*(4*x-1);
        N(:,10) = -(16*x.*(2*x-1/2).*(4*x-2))/3;
        N(:,5) = (16*x.*(2*x-1/2).*(4*x-2))/3;
        N(:,2) = 0;
    end
elseif strcmp(basis_type, "P5")
    lambda = 1-x-y;
    if strcmp(derivative_type, "x")
        N(:,1) = (5*lambda-4).*((5*lambda)/2-3/2).*((5*lambda)/3-2/3).*((5*lambda)/4-1/4).*lambda;
        N(:,15) = (25*y.*(5*lambda-3).*((5*lambda)/2-1).*((5*lambda)/3-1/3).*lambda)/4;
        N(:,12) = (25*y.*(5*lambda-2).*((5*lambda)/2-1/2).*(5*y-1).*lambda)/6;
        N(:,9) = (25*y.*(5*lambda-1).*(5*y-2).*((5*y)/2-1/2).*lambda)/6;
        N(:,6) = (25*y.*(5*y-3).*((5*y)/2-1).*((5*y)/3-1/3).*lambda)/4;
        N(:,3) = y.*(5*y-4).*((5*y)/2-3/2).*((5*y)/3-2/3).*((5*y)/4-1/4);
        N(:,4) = (25*x.*(5*lambda-3).*((5*lambda)/2-1).*((5*lambda)/3-1/3).*lambda)/4;
        N(:,16) = (125*x.*y.*(5*lambda-2).*((5*lambda)/2-1/2).*lambda)/3;
        N(:,21) = (125*x.*y.*(5*lambda-1).*(5*y-1).*lambda)/4;
        N(:,18) = (125*x.*y.*(5*y-2).*((5*y)/2-1/2).*lambda)/3;
        N(:,14) = (25*x.*y.*(5*y-3).*((5*y)/2-1).*((5*y)/3-1/3))/4;
        N(:,7) = (25*x.*(5*lambda-2).*((5*lambda)/2-1/2).*(5*x-1).*lambda)/6;
        N(:,19) = (125*x.*y.*(5*lambda-1).*(5*x-1).*lambda)/4;
        N(:,20) = (125*x.*y.*(5*x-1).*(5*y-1).*lambda)/4;
        N(:,11) = (25*x.*y.*(5*x-1).*(5*y-2).*((5*y)/2-1/2))/6;
        N(:,10) = (25*x.*(5*lambda-1).*(5*x-2).*((5*x)/2-1/2).*lambda)/6;
        N(:,17) = (125*x.*y.*(5*x-2).*((5*x)/2-1/2).*lambda)/3;
        N(:,8) = (25*x.*y.*(5*x-2).*((5*x)/2-1/2).*(5*y-1))/6;
        N(:,13) = (25*x.*(5*x-3).*((5*x)/2-1).*((5*x)/3-1/3).*lambda)/4;
        N(:,5) = (25*x.*y.*(5*x-3).*((5*x)/2-1).*((5*x)/3-1/3))/4;
        N(:,2) = x.*(5*x-4).*((5*x)/2-3/2).*((5*x)/3-2/3).*((5*x)/4-1/4);
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -(5*(5*lambda-4).*((5*lambda)/2-3/2).*((5*lambda)/3-2/3).*lambda)/4-(5*(5*lambda-4).*((5*lambda)/2-3/2).*((5*lambda)/4-1/4).*lambda)/3-(5*(5*lambda-4).*((5*lambda)/3-2/3).*((5*lambda)/4-1/4).*lambda)/2-5*((5*lambda)/2-3/2).*((5*lambda)/3-2/3).*((5*lambda)/4-1/4).*lambda-(5*lambda-4).*((5*lambda)/2-3/2).*((5*lambda)/3-2/3).*((5*lambda)/4-1/4);
        N(:,15) = -(25*y.*(5*lambda-3).*((5*lambda)/2-1).*((5*lambda)/3-1/3))/4-(125*y.*(5*lambda-3).*((5*lambda)/2-1).*lambda)/12-(125*y.*(5*lambda-3).*((5*lambda)/3-1/3).*lambda)/8-(125*y.*((5*lambda)/2-1).*((5*lambda)/3-1/3).*lambda)/4;
        N(:,12) = -(25*y.*(5*lambda-2).*((5*lambda)/2-1/2).*(5*y-1))/6-(125*y.*(5*lambda-2).*(5*y-1).*lambda)/12-(125*y.*((5*lambda)/2-1/2).*(5*y-1).*lambda)/6;
        N(:,9) = -(25*y.*(5*lambda-1).*(5*y-2).*((5*y)/2-1/2))/6-(125*y.*(5*y-2).*((5*y)/2-1/2).*lambda)/6;
        N(:,6) = -(25*y.*(5*y-3).*((5*y)/2-1).*((5*y)/3-1/3))/4;
        N(:,3) = 0;
        N(:,4) = (25*(5*lambda-3).*((5*lambda)/2-1).*((5*lambda)/3-1/3).*lambda)/4-(25*x.*(5*lambda-3).*((5*lambda)/2-1).*((5*lambda)/3-1/3))/4-(125*x.*(5*lambda-3).*((5*lambda)/2-1).*lambda)/12-(125*x.*(5*lambda-3).*((5*lambda)/3-1/3).*lambda)/8-(125*x.*((5*lambda)/2-1).*((5*lambda)/3-1/3).*lambda)/4;
        N(:,16) = (125*y.*(5*lambda-2).*((5*lambda)/2-1/2).*lambda)/3-(125*x.*y.*(5*lambda-2).*((5*lambda)/2-1/2))/3-(625*x.*y.*(5*lambda-2).*lambda)/6-(625*x.*y.*((5*lambda)/2-1/2).*lambda)/3;
        N(:,21) = (125*y.*(5*lambda-1).*(5*y-1).*lambda)/4-(125*x.*y.*(5*lambda-1).*(5*y-1))/4-(625*x.*y.*(5*y-1).*lambda)/4;
        N(:,18) = (125*y.*(5*y-2).*((5*y)/2-1/2).*lambda)/3-(125*x.*y.*(5*y-2).*((5*y)/2-1/2))/3;
        N(:,14) = (25*y.*(5*y-3).*((5*y)/2-1).*((5*y)/3-1/3))/4;
        N(:,7) = (25*(5*lambda-2).*((5*lambda)/2-1/2).*(5*x-1).*lambda)/6-(25*x.*(5*lambda-2).*((5*lambda)/2-1/2).*(5*x-1))/6+(125*x.*(5*lambda-2).*((5*lambda)/2-1/2).*lambda)/6-(125*x.*(5*lambda-2).*(5*x-1).*lambda)/12-(125*x.*((5*lambda)/2-1/2).*(5*x-1).*lambda)/6;
        N(:,19) = (125*y.*(5*lambda-1).*(5*x-1).*lambda)/4-(125*x.*y.*(5*lambda-1).*(5*x-1))/4+(625*x.*y.*(5*lambda-1).*lambda)/4-(625*x.*y.*(5*x-1).*lambda)/4;
        N(:,20) = (125*y.*(5*x-1).*(5*y-1).*lambda)/4-(125*x.*y.*(5*x-1).*(5*y-1))/4+(625*x.*y.*(5*y-1).*lambda)/4;
        N(:,11) = (25*y.*(5*x-1).*(5*y-2).*((5*y)/2-1/2))/6+(125*x.*y.*(5*y-2).*((5*y)/2-1/2))/6;
        N(:,10) = (25*(5*lambda-1).*(5*x-2).*((5*x)/2-1/2).*lambda)/6-(25*x.*(5*lambda-1).*(5*x-2).*((5*x)/2-1/2))/6+(125*x.*(5*lambda-1).*(5*x-2).*lambda)/12+(125*x.*(5*lambda-1).*((5*x)/2-1/2).*lambda)/6-(125*x.*(5*x-2).*((5*x)/2-1/2).*lambda)/6;
        N(:,17) = (125*y.*(5*x-2).*((5*x)/2-1/2).*lambda)/3-(125*x.*y.*(5*x-2).*((5*x)/2-1/2))/3+(625*x.*y.*(5*x-2).*lambda)/6+(625*x.*y.*((5*x)/2-1/2).*lambda)/3;
        N(:,8) = (25*y.*(5*x-2).*((5*x)/2-1/2).*(5*y-1))/6+(125*x.*y.*(5*x-2).*(5*y-1))/12+(125*x.*y.*((5*x)/2-1/2).*(5*y-1))/6;
        N(:,13) = (25*(5*x-3).*((5*x)/2-1).*((5*x)/3-1/3).*lambda)/4-(25*x.*(5*x-3).*((5*x)/2-1).*((5*x)/3-1/3))/4+(125*x.*(5*x-3).*((5*x)/2-1).*lambda)/12+(125*x.*(5*x-3).*((5*x)/3-1/3).*lambda)/8+(125*x.*((5*x)/2-1).*((5*x)/3-1/3).*lambda)/4;
        N(:,5) = (25*y.*(5*x-3).*((5*x)/2-1).*((5*x)/3-1/3))/4+(125*x.*y.*(5*x-3).*((5*x)/2-1))/12+(125*x.*y.*(5*x-3).*((5*x)/3-1/3))/8+(125*x.*y.*((5*x)/2-1).*((5*x)/3-1/3))/4;
        N(:,2) = (5*x-4).*((5*x)/2-3/2).*((5*x)/3-2/3).*((5*x)/4-1/4)+(5*x.*(5*x-4).*((5*x)/2-3/2).*((5*x)/3-2/3))/4+(5*x.*(5*x-4).*((5*x)/2-3/2).*((5*x)/4-1/4))/3+(5*x.*(5*x-4).*((5*x)/3-2/3).*((5*x)/4-1/4))/2+5*x.*((5*x)/2-3/2).*((5*x)/3-2/3).*((5*x)/4-1/4);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -(5*(5*lambda-4).*((5*lambda)/2-3/2).*((5*lambda)/3-2/3).*lambda)/4-(5*(5*lambda-4).*((5*lambda)/2-3/2).*((5*lambda)/4-1/4).*lambda)/3-(5*(5*lambda-4).*((5*lambda)/3-2/3).*((5*lambda)/4-1/4).*lambda)/2-5*((5*lambda)/2-3/2).*((5*lambda)/3-2/3).*((5*lambda)/4-1/4).*lambda-(5*lambda-4).*((5*lambda)/2-3/2).*((5*lambda)/3-2/3).*((5*lambda)/4-1/4);
        N(:,15) = (25*(5*lambda-3).*((5*lambda)/2-1).*((5*lambda)/3-1/3).*lambda)/4-(25*y.*(5*lambda-3).*((5*lambda)/2-1).*((5*lambda)/3-1/3))/4-(125*y.*(5*lambda-3).*((5*lambda)/2-1).*lambda)/12-(125*y.*(5*lambda-3).*((5*lambda)/3-1/3).*lambda)/8-(125*y.*((5*lambda)/2-1).*((5*lambda)/3-1/3).*lambda)/4;
        N(:,12) = (25*(5*lambda-2).*((5*lambda)/2-1/2).*(5*y-1).*lambda)/6-(25*y.*(5*lambda-2).*((5*lambda)/2-1/2).*(5*y-1))/6+(125*y.*(5*lambda-2).*((5*lambda)/2-1/2).*lambda)/6-(125*y.*(5*lambda-2).*(5*y-1).*lambda)/12-(125*y.*((5*lambda)/2-1/2).*(5*y-1).*lambda)/6;
        N(:,9) = (25*(5*lambda-1).*(5*y-2).*((5*y)/2-1/2).*lambda)/6-(25*y.*(5*lambda-1).*(5*y-2).*((5*y)/2-1/2))/6+(125*y.*(5*lambda-1).*(5*y-2).*lambda)/12+(125*y.*(5*lambda-1).*((5*y)/2-1/2).*lambda)/6-(125*y.*(5*y-2).*((5*y)/2-1/2).*lambda)/6;
        N(:,6) = (25*(5*y-3).*((5*y)/2-1).*((5*y)/3-1/3).*lambda)/4-(25*y.*(5*y-3).*((5*y)/2-1).*((5*y)/3-1/3))/4+(125*y.*(5*y-3).*((5*y)/2-1).*lambda)/12+(125*y.*(5*y-3).*((5*y)/3-1/3).*lambda)/8+(125*y.*((5*y)/2-1).*((5*y)/3-1/3).*lambda)/4;
        N(:,3) = (5*y-4).*((5*y)/2-3/2).*((5*y)/3-2/3).*((5*y)/4-1/4)+(5*y.*(5*y-4).*((5*y)/2-3/2).*((5*y)/3-2/3))/4+(5*y.*(5*y-4).*((5*y)/2-3/2).*((5*y)/4-1/4))/3+(5*y.*(5*y-4).*((5*y)/3-2/3).*((5*y)/4-1/4))/2+5*y.*((5*y)/2-3/2).*((5*y)/3-2/3).*((5*y)/4-1/4);
        N(:,4) = -(25*x.*(5*lambda-3).*((5*lambda)/2-1).*((5*lambda)/3-1/3))/4-(125*x.*(5*lambda-3).*((5*lambda)/2-1).*lambda)/12-(125*x.*(5*lambda-3).*((5*lambda)/3-1/3).*lambda)/8-(125*x.*((5*lambda)/2-1).*((5*lambda)/3-1/3).*lambda)/4;
        N(:,16) = (125*x.*(5*lambda-2).*((5*lambda)/2-1/2).*lambda)/3-(125*x.*y.*(5*lambda-2).*((5*lambda)/2-1/2))/3-(625*x.*y.*(5*lambda-2).*lambda)/6-(625*x.*y.*((5*lambda)/2-1/2).*lambda)/3;
        N(:,21) = (125*x.*(5*lambda-1).*(5*y-1).*lambda)/4-(125*x.*y.*(5*lambda-1).*(5*y-1))/4+(625*x.*y.*(5*lambda-1).*lambda)/4-(625*x.*y.*(5*y-1).*lambda)/4;
        N(:,18) = (125*x.*(5*y-2).*((5*y)/2-1/2).*lambda)/3-(125*x.*y.*(5*y-2).*((5*y)/2-1/2))/3+(625*x.*y.*(5*y-2).*lambda)/6+(625*x.*y.*((5*y)/2-1/2).*lambda)/3;
        N(:,14) = (25*x.*(5*y-3).*((5*y)/2-1).*((5*y)/3-1/3))/4+(125*x.*y.*(5*y-3).*((5*y)/2-1))/12+(125*x.*y.*(5*y-3).*((5*y)/3-1/3))/8+(125*x.*y.*((5*y)/2-1).*((5*y)/3-1/3))/4;
        N(:,7) = -(25*x.*(5*lambda-2).*((5*lambda)/2-1/2).*(5*x-1))/6-(125*x.*(5*lambda-2).*(5*x-1).*lambda)/12-(125*x.*((5*lambda)/2-1/2).*(5*x-1).*lambda)/6;
        N(:,19) = (125*x.*(5*lambda-1).*(5*x-1).*lambda)/4-(125*x.*y.*(5*lambda-1).*(5*x-1))/4-(625*x.*y.*(5*x-1).*lambda)/4;
        N(:,20) = (125*x.*(5*x-1).*(5*y-1).*lambda)/4-(125*x.*y.*(5*x-1).*(5*y-1))/4+(625*x.*y.*(5*x-1).*lambda)/4;
        N(:,11) = (25*x.*(5*x-1).*(5*y-2).*((5*y)/2-1/2))/6+(125*x.*y.*(5*x-1).*(5*y-2))/12+(125*x.*y.*(5*x-1).*((5*y)/2-1/2))/6;
        N(:,10) = -(25*x.*(5*lambda-1).*(5*x-2).*((5*x)/2-1/2))/6-(125*x.*(5*x-2).*((5*x)/2-1/2).*lambda)/6;
        N(:,17) = (125*x.*(5*x-2).*((5*x)/2-1/2).*lambda)/3-(125*x.*y.*(5*x-2).*((5*x)/2-1/2))/3;
        N(:,8) = (25*x.*(5*x-2).*((5*x)/2-1/2).*(5*y-1))/6+(125*x.*y.*(5*x-2).*((5*x)/2-1/2))/6;
        N(:,13) = -(25*x.*(5*x-3).*((5*x)/2-1).*((5*x)/3-1/3))/4;
        N(:,5) = (25*x.*(5*x-3).*((5*x)/2-1).*((5*x)/3-1/3))/4;
        N(:,2) = 0;
    end
end
end